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Abstract 



In 1917, Johann Radon introduced the Radon transform [33] which is used in 1963 by Allan 
MacLeod Cormack for application in the context of tomographic image reconstruction. He pro- 
posed to reconstruct the spatial variation of the material density of the body from X-Ray images 
(radiographies) for different directions; He implemented this method and made a test for a proto- 
type Computarized Tomography (CT) scanner [3]. Independently, Godfrey Newbold Hounsfield 
derived an algorithm and built the first medical CT scanner. This was a great achievement for 
the twentieth century, because one can see inside an object without opening it up; Cormack and 
Hounsfield won the Nobel Prize of Medicine in 1979 for the development of computer assisted 
tomography. 

Basically the idea of the X-ray CT is to get an image of the interior structure of an object by 
X-raying the object from many different directions. X rays go in straight lines inside the body and 
its energy is attenuated more and less depending on the density of the matter in his trajectory. So 
the simplest model relating the log ratio (ln///o) of the observed energy / with emitted energy 
Iq to the spatial spatial distribution / of the body is a line integral. The mathematical problem 
is then estimating a multivariate function / from its line integrals. 

Four year before Cormack's idea, Abe Sklar introduced another theory in the context of Statistics 
called «copula». He gave the theorem which now bears his name [35j. Succinctly stated, copulas 
are functions that link multivariate distributions to theirs univariate marginal functions. 

One of the problems arising in Statistics is the reconstruction of joint distribution function from its 
given marginal functions. It appeared that copulas captivated all dependence structure concerning 
the marginal functions and offer a wide range of parametric family model which could be used as 
a model for the joint distribution function. This problem is the same as in Tomography, because 
a marginal density is obtained from a line integral of its joint distribution. The analogy is then 
that the joint density / has to be reconstructed from its marginals who are obtained by their 
integration over lines. 

To achieve our goal based on the mathematical aspects of tomography in imaging sciences using 
copulas, we give some prerequisites about copulas and tomography. In the particular case of only 
given horizontal and vertical projections corresponding to a given two marginal functions, we link 
the theory of copula to tomography via the Radon transform and Sklar's theorem. Let us note 
that determining the density functions (or the object) from two projections is an ill-posed inverse 
problem. 

Finally, to simulate an image reconstruction we use some members of copulas families (Archimedean, 
Elliptic) already in widespread use. We have also written a package so called «copula-tomography» 
to handle all those copulas and allow users to simulate a tomographic image reconstruction 
through copulas. Some preliminary results for a few number of projections are promising. 
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Introduction 



The word «copula» originates from the Latin meaning «link, chain; union». In statistical liter- 
ature, according to the seminal result in the copula's theory stated by Abe Sklar [35j in 1959; 
a copula is a function that connects a multivariate distribution function to its given univariate 
marginal distributions. 

The main point of our report, is bridging the gap between the theory of copula and tomography. 
We consider the inverse problem which is the reconstruction in Computed Tomography (CT) 
of an object using the Radon transform from only two horizontal and vertical projections. This 
problem of reconstruction is similar, in the statistical point of view, when it is known or assumed 
to know the distributions of each random variables but not their joint distribution function, or 
their copula density. 

There is an increasing interest concerning copulas, widely used in Financial Mathematics [23], in 
modelling of Environmental Data [20j. Recently, in Computational Biology, copulas are used for 
the reconstruction of accurate cellular networks [18j. Copula appeared to be a new powerful tool 
to model the structure of dependence. Copulas are useful for constructing joint distributions, 
particularly with nonnormal random variables. 

What are copulas and how they are related and suitable in the modelling processes of image 
reconstruction? 

To give more details in order to answer those questions, we organise our report as follows: in 
the next chapter, we recall some definitions related to multivariate distribution functions and we 
give the Sklar's theorem, highlighted by some methods to generate a new copula. In the chapter 
3, through many illustrations, we show some parametrics family of copulas. We dedicated the 
chapter 4 on operation yielding to discrete copula via bistochastic matrices. We expose some well 
known statistical methods such that the Maximum Likelihood Estimation (MLE), the Inference 
Functions for Margins (IFM) and the Bayesian estimation, in the chapter 5. And we give also 
other methods from literature dealing with the good way to choose the right copulas. In the 
chapter 6, we start by explanation of the basic mathematical model of X-ray tomography. The 
top point is discussed in the chapter 7, where we link the theory of copula to tomographic image 
reconstruction through the Radon Transform. We focus on the case of only two given horizontal 
and vertical projections which correspond to a given two marginal functions. This is an ill-posed 
inverse problem because information about the density functions (or the object) to reconstruct is 
obtained from indirect and limited data. 

Andrei Nikolaevich Tikhonov (Russian mathematician, 1906-1993) said: 



1 



Page 2 



For a long time mathematicians felt that ill-posed problems cannot describe 
real phenomena and objects. However [...] the class of ill-posed problems 
includes many classical mathematical problems and, most significantly, 
that such problems have important applications. 

For our goal of image reconstruction, some Archimedean family of copula, the Elliptic class of 
copula in particular the Gaussian copula and mixture of Gaussian density, are used successively for 
simulation. We present some results we have obtained and a summary of our future work in the 
last chapter, also a short description concerning the «copula-tomography» package in appendix. 
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Copula in Statistics 



2.1 Multivariate distribution 

In this section, we first give a short summary of definition and properties of a multivariate 
distribution. We will then focus on copula. We extend the notion of increasing function of one 
variable to n variables, so called n-increasing function, then we introduce the Sklar's theorem 
and also a method to generate a copula. We assume that our random variables are continuous if 
necessary. 

2.1.1 Joint probability density function(pdf) 

Let Xi, . . . , Xnhe n continuous random variables0all defined on the same probability space. The 
joint probability density function (pdf) of denoted by /(Xi, ,Xn)(^i5 ^^n). 

is the function / : ^ M such that for any domain C in the n-dimensional space of 
the values of the variables Xi, . . . , X^, the probability that a realisation of the set variables falls 
inside the domain D is 

Prob(Xi,X2,...,X^ G i^) = / f{xi,...,xn)dxi...dxn. (2.1) 

JD 

Property 2.1. The two following properties are satisfied: 

• / f{Xi,...,Xn)dXi ...dXn = l. 

2.1.2 Cumulative distribution function(cdf) 

The cumulative distribution function (cdf) is defined by the formula : 

... / f{yi,...,yn)dyi ... dyn- 

-oo J — oo 

■"^Conciseiy a n-dimensional random variables is a function from a sample space S of an experiment into M.^. 
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There is a relation between the pdf and the cdf: 

2.1.3 Gaussian Mixture(GM) distribution 

The Gaussian mixture (GM) model is one of the most used model in statistics and modelling 
process (for clustering, or density estimation). It is defined via a Gaussian (or normal) distribution. 

Let X = be a n-dimensional random vector that is multivariate normally dis- 

tributed, then the probability density function (pdf) is 

^^"^ I ^' ^ (27r)i|S|V2 {-I - - ^)) ' (2.3) 

where 

• fi — [/ii, . . . , /in]^, is the vector of the mean values /i^ = E(X^), 

• S is the covariance matrix, a n x n non-singular, positive definite real matrix, with entries 

j:,j^E[{Xi-f,i){Xj-ixj)]. 

When the distribution of X is the multivariate normal distribution, we will use the following 
notation: 



One could also derive the cumulative distribution function (cdf), which is simply 

"J (27r)"/^ I S |V2 ^^P[~2 (^ ~ ^"^ (^ ~ • • • 



(2.4) 



In the 2-dimensional case, the pdf has the form: 

/(^i' ^2) = 7== exp —- — + 



where p is the correlation between xi and X2, and S = 



For the GM distribution, each of the K components X^, is such that X^ ~ J^ifik, ^k), with the 
probabilistic component weights ak- Therefore, the GM pdf has the following form in the case 
of discrete variables (change the sum to an integral for continuous variables): 

K 

/((c) = ^afcAr(/Xfc,Sfc), (2.5) 



k=i 

K 



where < < 1 , aj. = 1, and = (xi, . . . , x^). 
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2.1.4 t-distribution 

The t— distribution (or the Student's t-distribution) with v degrees of freedom has the probability 
density function : 

= ^7^(1 + (2-6) 



^B(|,i)V^ ■ y 

where B(a,/?) is the beta function which is defined as. 

B{p,q) = [ xP-\l-xy-^dx 
Jo 

for any real numbers p^q > 0. We will denote X ^ t{iy). We could also write the beta function 
in term of the gamma function, 

poo 

T{z) = / e-H'-^ dt 



for z G C with > 0, and by analytic continuation for the rest of the complex plane, except 
for the non-positive integers, where it has simple poles. Another equivalent definition is 



n=l 



where 7 is Euler's constant. 

Property 2.2. Remind about properties of the beta and gamma functions: 

B{p, q) = ^^J^^^^'^^ for all complex numbers p and q for which the right-hand side is defined. 

V{p + q) 

AlsoBip,q) = B{q,p) and B (l,^] = w. 



r 2^ 

T{p) = {p — l)r(p — 1) and for any integer n> 1, T{n) = (n — 1)! 
See [1^ for more details. 

Therefore, the probability density function of the Student's t-distribution is : 

/(^) = ^=^ 1 + - ■ (2.7) 



^r(i) V y 

Property 2.3. l/l/e recall some useful properties of the t-distribution: 

• t(l) = Cauchy(0, 1), the Cauchy distribution with parameters and 1. 

• t{iy) ^ A/'(0, 1) when v — > oo. 



IfX - t(zy), E[\X 1^] exists if and only if k < v. 
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The last property means that, there is no mean for a t— distribution \f u = 1. And for k = 1, that 
is z/ > 1, E[X] = 0; then for the case u > 2, the variance of X is equal to If Xi, . . . , 
are random samples from a normal distribution with mean /i and variance a^. If we denote JI, 

n n 

TT _ 1 ^ ^^A 4.u^ u ^2 ^ 1 -\2 

k=l k=l 



respectively the sample mean M = ^ ^"^^ ^'^^ sample variance = ^ (x^ — /x) , then 



jji — jji 



(J 



(2.8) 



One can define the multivariate t-distribution. Let X be a vector of n variate t-distribution with 
V degrees of freedom and mean vector /x, according to the previous property in the univariate 
case, we set the mean /x > 1 and the degree of freedom v > 2, therefore the covariance matrix 



is ^S. 

v—2 



Hence X takes the following form 



X 



Vs 



where Z ^ A/'(0, E) and S is distributed independently of Z and has a xl distribution, then the 
multivariate t-distribution is given by 



-1/2 



l + ^{x- At)^S-^ {x - fj,) 



-(^) 



(2.9) 



From the graphs below with u = 2, u = 20, we could clearly see that as u increases, asymptotically 
the t— distribution reaches the standard normal distribution. 





2.1.5 n-dimensional marginal distribution 

Definition 2.4. Let {X^ : i G /} be a set of random variables such that the subset I C 
{1, 2, . . . , n} . If we denote F (xi, . . . , Xn) the joint cdf defined on R^. The marginal distri- 
bution F^Xi-.iei}' obtained by summing (for discrete variables) 

or integrating (for continuous variables) 

F{X,:ieI} (Xi) = Fi (Xi) = / • • • , ^n) Yl 

over all values of the other variables. 
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Remark Fi{xi) = F(xi,oc) and F2{x2) = F(oc,X2), for continuous bivariate distribution. 




2.2 Sklar's Theorem 



2.2.1 Multivariate Copula 

A n— copula (or copula) is a multivariate joint distribution defined on the n-dimensional unit cube 
[0, 1]^ such that every marginal distribution is uniform on the interval [0, 1] . 

Definition 2.5. A multivariate copula, or a n— copula denoted by C , is a function from [0, 1]^ 
to [0, 1] with the following properties : 

1. Vu=K,...,^OG[0,ir 

• C{ui, . . . , Un) = 0, if at least one component of u is equal to zero , 
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2. Vu=(Mi,...,M„)G[0,ir, 

• C(l, . . . , 1, Mj, 1, . . . , 1) = Mj, if all components of u are equal to 1 except Ui ; 

3. V(zii'\ . . . , u^n ^) e [0, 1]", anc/ V(zif \ . . . , e [0, 1]" sue/? t/?at tif ^ < «f ^ Vj 



2 ^ ^ 



2 

Zl = l Zn = l 



7"/7e last property below is a n dimensional analogue version of an univariate increasing function. 
Therefore any copula . . . , Un) is n— increasing function. We can also express this idea in 

the following way. 



2.2.2 C-Volume 

Let a = (ai, a2, . . . , a^) and 6 = (61, 62, • • • , &n)- An n-Box, denoted by B is the Cartesian 
product [ai, hi] x [a2, 62] x • • • [(^td &n]- ^Iso denote by c = (ci, C2, . . . , c^) the vertices of B, 
i.e. Ck is equal to ak or hk, and a < b be equivalent to ak < bk V fc. 

Definition 2.6. Let Si, S2, • • • , Sn be nonempty subsets ofR, and C an n— place real 
function^ such that Dom C = Si x S2, . . . , Sn- Let B [a, b] be an n—box such that all of its 
vertices are in DomC. Then the C— volume of B is given by 

Vc{B) = J2s9n{c)C{c), (2.10) 

ceB 

where 

sgn{c) = 



1, if Ck ak for an even number of k's 
— 1, if Ck = ak for an odd number of k's. 



We can also express, the (7— volume of an n—box B in the term of the nth order difference of 
C ox\ B 

Vc{B) = A3^C(t) = A'z^t--\ • • • ^'al<C{t), 
where the n first order difference of C is 

^afc^(t) = C {ti, . . . , tk-l, bk, • • • ,tn) — C {ti, . . . , tk-l, ak, • • • ,tn) • 

Definition 2.7. C is an n-increasing function if the C— volume of B, Vc{B) > 0. 

A copula C induces a probability measure on [0, 1]^ via Vc'([0, 1] x [0, 1] . . . [0, 1]) = C{ui, . . . , Un) 
(Caratheodory's theorem (measure theory)) see [34j. As consequence of definition I2.7[ any copula 
C satisfies the following inequality (see [30] ). 



^function whose domain, DomC, is a subset of R"^ and whose range, RanC, is a subset of '. 
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Theorem 2.8. Let C be an n-copula. Then for every u = (t^i, . . . , Un) and v = {vi, . . . ,Vn) in 

[0, ir , 

n 

\C{u)-Civ)\<^\uk-Vk\. 

k=0 

Therefore C is uniformly continuous on [0, 1]^ . 



2.2.3 Bivariate Copula 

A bivariate copula, or shortly a copula is a function from [0, 1]^ to [0, 1] 
Property 2.9. with the following properties from the previous with n = 2. 

1. V^, V G [0, 1] , C{u, 0) = = C(0, v); 

2. V^, ^ G [0, 1] , C{u,l) = u and C{l,v)=v; 

3. \/ui^U2^Vi^V2 G [0, 1] such that Ui < U2 and Vi < 

Vc {[Ui, U2] X [^1, ^2]) = C{U2, V2) - C{U2, vi) - C{ui, V2) + C{ui, Vi) > 0. 

To show the last property in the bivariate case, one have to set u^"^ = Ui and 2/2^ = i ^ 
{1, 2} in the last property for the above multivariate case. And to make sure that this relation is 
the generalisation of the definition of an increasing univariate function, one can set for example 

^1 = and V2 1. 

2.2.4 Sklar's Theorem 

The Sklar's theorem watertight the theory of copula. The proof of the 2-dimensional case can 
be found in [35]. We refer to [36j, for the proof of the n— dimensional version of this theorem. 
It is the most important result concerning copulas and widely used in all applications. 

Theorem 2.10. (Sklar's Theorem) 

Let F be a two-dimensional cumulative distribution function with marginal distributions functions 
Fi and F2. Then there exists a copula C such that: 

F(xi,X2) = C(Fi(xi),F2(x2)). (2.11) 

Furthermore, if the marginal functions are continuous, then the copula C is unique, and is given 
by 

CK,M2) = F{F{\m),F^\u2)). (2.12) 

Otherwise, if Fi or F2 is discontinuous, then C is uniquely determined on RanFi x RanF2. 
Conversely, for any univariate distribution functions Fi and F2 and any copula C, the function 
F is a two-dimensional distribution function with marginals Fi and F2, given by 1^2. 
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Along this report, relation (I2.12p appeared to be very useful in practice. 
Theorem 2.11. (n-dimensional Sklar's Theorem) 

Let F be a joint distribution function with marginals cumulative distribution functions Fi, . . . , F^. 
Then there exists a copula C such that for all xi, . . . ,Xn 

F(xi, . . . , = . . . , Fd{xn)). (2.13) 

Furthermore, //Fi, . . . , are continuous functions, then the copula C is unique and for all 

CK, . . . , ^0 = F{F,-\u,), . . . , F-^K)). (2.14) 

Otherwise, if there is at least one marginal discontinuous, then C is uniquely determined on 
RanFi x • • • x RanF^. Conversely, suppose C is a copula and that Fi{ui), . . . , Fn{un) are 
the univariate cumulative distribution functions. The function F defined as follows: 

F(xi, . . . , = C(Fi(xi), . . . , Fa{xn)) (2.15) 

is a joint distribution function with marginals cumulative distribution functions Fi{ui),. . . ,Fn{un)- 

From the Sklar's Theorem, one can understand already why copula can be seen as a powerfull 
tools in modelling dependences of several random variables. 



2.2.5 Copula Density 

From (12. 2p and differentiating (I2.14P gives the density of a copula 

d-C f[F,-\u^),...,F-\u^)] 



dui . . . dur, 



(2.16) 



i=l 



where / is the joint probability density function of the cdf F and each fi is the marginal density 
functions of the marginal cdf F^. 
The pdf could be expressed as : 

n 
i=l 

Equivalently, since Xj = F-~^(Mj) — > Ui = Fi{xi), it follows that 



f{Xi,...,Xn)=C {Fi{Xi), Fn{Xn)) JJ fi (Xi) 



(2.17) 



Equation (I2.17p is central to attend the goal of this report. 



Section 2.3. Generating Copulas by the Inversion Method 



Page 11 



2.3 Generating Copulas by the Inversion Method 

Given Xi and X2 two random variables. F(xi,X2) a joint distribution function and its margins 
Fi(xi) and F2{x2), all assumed to be continuous. The corresponding copula can be constructed 
using the unique inverse transformations (Quantile transform) Xi = F^^{Ui), X2 = F2^{U2), 
where Ui, and U2 are uniformly distributed on [0, 1]. 

This method based directly on Sklar's theorem, generates a copula from the following equation 

C{ui,U2) = F {F{^ {ui) , {U2)) = F(xi,X2), where i^i, 2/2 ai'e uniform on [0,1]. (2.18) 

Now let us give an example to illustrate, how to construct a parametric family of copula using 
inversion method. 

Starting with the joint distribution function: 

F^(xi, X2) = [(1 + e-^^ + + (1 - a)e-^^-^^] for Xi, X2 G 1. 

By taking limit when X2 and xi goes to infinity successively in the above expression leads to the 
marginal distribution functions, which are 

Fi(xi) = F,(xi,oc) = (1 + 6-^^-1 

and 

F2{x2) = F^{^,X2) = {l + e-^-)-\ 



Simple algebra gives us the inverse transformations 



and 



X2 = F^\u2) = \n 



Finally after substituting in (I2.18p yields to: 

1 — Ui I — U2 



C{ui,U2) 
C{ui,U2) 

C{ui, U2) 



1 + 



Ui 



+ 



U2 



+ {l-a) 



{l-Ur){l-U2) 



U1U2 



1 — a{l — Ui){l — U2) 



U1U2 
U1U2 



1 — a{l — Ui){l — U2) 



where a G [—1, 1) . 



(2.19) 



The copula (I2.19p we have constructed is known as the Ali-Mikhail-Haq copula. It was introduced 
in 1978, see Ali et a/.[lj. 
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Remark Construction of new families of copulas by different other existing methods (algebraic 
or geometric) is an important domain of research. For example as a tool of constructing a new 
copula, the geometric method takes account only of the definition of copula, without using any 
given distribution function. In order that the third property (see Definition 12. 7p about copula 
holds, prior knowledge about the geometric representation is used. Those informations could be 
the support of the function, the graphical shape of its different sections (horizontal, vertical or 
diagonal). 
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Some Families of Copulas 



There are many copulas and each of them has a specific property. In this section, we list few of 
the widely used families. We will end by definition of discrete copula associated to a bistochastic 
matrix, and some statistical methods followed by discussion on the way to choose the right copula. 
Some graphical representations of copulas are listed in Appendix [Bl 



3.1 Usual Copulas 

The product copula (or independent copula) is the simplest copula, has the form 

n(t^i, . . . , Un) = Ui . . .Un for all Ui G [0, 1] , 

corresponds to independence, therefore it is important as a benchmark. 

The Frechet-Hoeffding upper bound copula (or comonotonicity copula) is : 

M(t^i, . . . , Un) = min {t^i, . . . , Un} for all Ui G [0, 1] . 
And the Frechet-Hoeffding lower bound (or countermonotonicity copula) : 

W{ui, . . . , Un) = max < 1 — n + ^ Ui, > for all Ui G [0, 1] . 



i=l 



One can easily check that all properties about copula are satisfied by W, U and M. 

Property 3.1. Any copula C{ui, . . . ,Un), satisfied the inequality called the Frechet-Hoeffding 
bound inequality 

W{Ui, ...,Un)< C{Ui, ...,Un)< M{ui, . . . , ^n)- 

Remark The upper bound is a copula for any value of n and the lower bound is a copula for 
n = 2. For n > 2 the lower bound may be a copula under some condition given in Theorem 3.6 
from the book t21j. 
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3.2 Archimedean Copulas 



The Archimedean copulas is one of the important classes of copulas which has many applications. 
More details about this section can be found in [30j page 89. 

Definition 3.2. Let Lp be a continuous, strictly decreasing function from [0, 1] to [0, oc] such 
that Lp{l) 0. The pseudo-inverse of Lp is the function Lp^~^^ with DomLp^~^^ [0, oc] and 

RanLp^~^^ = [0, 1] given by 



(p-^{t) 0<t<(^(0), 
^(0) < t < oo. 



Note that ip^ ^1 is continuous and decreasing on [0, oo], and strictly decreasing on [0, (/^(O)] 
Furthermore, ip^~'^^ {ip{ui)) = ui on [0, 1] , and 



t 0<t< (p{0), 

ip{0) (p{0) <t<oo. 

min(t,(/7(0)) . (3.1) 



Finally, if(p{0) = oo, then ip^ = ip ^. 



Theorem 3.3. Let ip be a continuous, strictly decreasing function from [0, 1] to [0, oo] such that 
(p{l) = 0, and let (p^~^^ be the pseudo-inverse of tp. Let C be the function from [0, 1]^ to [0, 1] 
given by 

C(tii, «2) = ^'-'1 {.^{ui) + ip{u2)) . (3.2) 
Then C is a copula if and only if Lp is convex. 
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Archimedean copulas are in the form (13. 2p and the function is called the generator of the 
copula. 99 is a strict generator if 99(0) = oc, then 99'"^! = Lp~^ and the relation 

gives a strict Archimedean copula. 

Property 3.4. The following algebraic properties are satisfied by any Archimedean copula C , 
those properties distinguish this class of copula from all other copula. 

1. C{ui,U2) C{u2, ui) meaning that C is symmetric M ui,U2 G [0, 1] ; 

2. C is associative W ui^U2^ Us G [0,1] i.e. C{C{ui^U2)^us) = C{ui^C{u2^us))] 

3. lfa>0 is any constant then a^p is again a generator of C . 

From the relation (13. 2p . we clearly see that all information about the dependence structure of the 
multivariate copula is reduced only to the study of the univariate generator Lp. This characteristic 
of the Archimedean family made them more attractive. 



The following theorem gives a technique to find a generator of an Archimedean copula. 



Theorem 3.5. Let C be an Archimedean copula with generator (f in Q. Then for almost all ui 

and U2 in [0, 1], 

/ dC{ui,U2) /. dC{ui,U2) f^^. 



Definition 3.6. If F{xi,X2, • • • ,Xn), and Fi{xi) denoted respectively the multivariate distribu- 
tion and its marginal functions, one particularly simple form of a n— dimensional Archimedean 
is 



where Lp is the generator function such that Lp{l) 0, Lp{Q) oc; and satisfies the convexity 
properties 99 (x) < 0, 99 (x) > 0. 

One easy way to compute the bivariate copula density function c{ui^U2) of the copula C{ui^U2), 
using the generator function (p under some conditions is given by: 

[if'{C{ui,U2))] 

Remark Other rigorous mathematics way to define the Archimedean copula is related to the 
Laplace transform (for details and beauty of this method, we refer to [27j). 
Let A be a distribution function with support IR+ and (p its Laplace transform, 



Jo 



exp(— tx) A(rfx), 
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Lp is strictly nondecreasing function, 99(0) = 1, Lp{-\-oo) = 0, then the following relation define a 
copula 

3.2.1 Independent copula 

For all ui^ U2 G [0, 1] , n(2ii, U2) = uiU2^ with the generator Lp{x) = — ln{x). And its copula 
density is identically equal to one. 

3.2.2 Gumbel copula 

This copula was originally studied by Gumbel in 1960 ( see [T4]) 

Ca{ui,U2) = exp [(-Int^i)"^ + {-lnu2)'^]^^ , 

with the generator Lp{x) = (— In(x))'^ and a G [1, oc) . 
And the density has the following form: 

Ca{Ui,U2) = 

(-ln2ii)-^+" [-l + a+(-ln2ii)" + (-ln2i2)"]^ [(- In^ii)" + (- In2i2)1"'+^ (- In^is)"'^" 
1 

X 



exp{(-ln2ii)^ + {-lnu2)'^y^'^ U1U2 

3.2.3 Clayton copula 

The following copula family was discussed in 1978 by Clayton [5j : 

Ca{ui,U2) = max (^[^i^"" + u^"" - l]^ ,0^ , a e [-1, oc)\{0} . 

The generator is Lp{t) = ^ {t~^ — 1). 
We compute its density which is given by: 

Ca{Ui,U2) = {l + a) 2/^'""^^'"" (-1 + + ^r)"'"^ • 

For a = the random variables are statistical independent, since Co{ui,U2) = H. 

But for large value of a, we have Coo(^i,^2) = M and the limiting case a = —1, yields to 

C{ui,U2) = W. 
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3.2.4 Ali-Mikhail-Haq copula 

Ali-Mikhail-Haq family: 

and the generator function is (f{x) = ln( "^~^^"^~^^ ). Its density is given by : 

-1 + a'^ (—1 + Ui + U2 — U1U2) — a (—2 + Ui + U2 + U1U2) 



Ca{Ui,U2) = 



[-1+ a (-1+^1) (-1 + ^2)]' 



3.2.5 Frank copula 

Frank copula has the following form (see [9]): 

Caiui,U2) = -^In (yZ^ - ^"") - ^1 - - e-«"^)]) , 

where a ^ and the generator is Lp{x) = In { ^^Iz^ ) • 
Its density is 

aexp [a (1 + t^i + 1^2)] [-1 + exp(a)] 

^a\^l)^2) — 2' 

{exp(a) — exp [a{l + 2/1)] — exp [a{l + 2/2)] + exp [a{ui + 1^2)]} 

3.2.6 Joe Copula 

Joe copula [19j has the form of 

Ca{ui,U2) = 1 - [(1 - ^1)" + (1 - U2T - (1 - ^1)" (1 - ^2)1^ , where a G [1, oc) , 
and the generator function is ip{x) = — ln(l — (1 — x)^). 
Its corresponding density is : 

Ca(tii, U2) = (1 - {a - [-1 + (1 - u,r] [-1 + (1 - u^T]} 

X [(1 - + (1 - «2)" - (1 - Ml)" (1 - «2)1"'+^ (1 - «2)"'+" . 

3.3 The Farlie-Gumbel-Morgenstern family 

Farlie-Gumbel-Morgenstern (shortly "FGM") copula was introduced in the basic functional follow- 
ing form by Eyraud [8J in 1938 and it was also discussed by Morgenstern ^22] in 1956 : 

Ca{ui, U2) = U1U2 (1 -l- a{l — ui){l — U2)) , where — 1 < a < 1. 
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For a 0, the FGM copula felt to be associative, then this family is not an Archimedean copula. 
One can check, for example 

FGM copula could be seen as a perturbation of the product copula, because the special case 
a = 0, leads to 11 the only Archimedean member. 

Its density is given by 

Ca{ui,U2) = 1 + a{l - 2ui){l - 2U2). 

3.4 Copula with cubic section 

Copula with cubic section is defined by Nelsen [30j and has the form: 

CaA^i^ ^2) = U1U2 + uiU2{l - ui){l - U2) [(3 + a{l - 2ui){l - 22x2)] , 
where are real constants such that a G [—1,2] and satisfy the conditions 
I /? 1^ a + 1 for a G [-1, 1/2] , and | /? |^ (6a - ?>a^)^ for a G [1/2, 2] . 
In order to have a nice representation, we consider the simple case where /5 = 0, 

Ca{ui, U2) = U1U2 + auiU2{l - ui){l - U2){1 - 2ui){l - 2U2). 
Its density is given by 

Ca{ui^U2) = 1 + a{l — 6ui + 6ul){l — 6u2 + 6ul). 

3.5 Elliptical copulas 

Elliptical copulas are the copula of elliptical distributions. We consider the examples of Gaussian 
distribution and t-distribution. 

3.5.1 Gaussian Copula 

The multivariate Gaussian copula is the copula of n-dimensional random vector that is 
multivariate normally distributed. This copula was proposed in 1983 by Lee t25j 

. . . , ^0 = *E • • • , *"'K)) , (3.5) 

^meaning that for given ui the function C{ui,U2) is cubic in ui and vice versa. 
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where $5] is the multivariate cdf (12. 4p and 

H^) = j exp ) dt (3.6) 

The bivariate Gaussian copula with the covariance matrix. 

-(;;: 

with Pearson's product-moment correlation coefficient p, therefore we have explicitly 

The limits cases p = —1, 0, 1 correspond respectively to the W{ui, U2) , Il{ui, U2) and M{ui, U2) 
copulas. Let us show the easy case p = 0. 

Cp{ui, U2) = / — exp <^ ^ dx dy 



dy 



In the multivariate case, from (I2.16P the multivariate Gaussian copula density is 

c,{ui, . . . , ^0 = |Sr^ exp 1^*^ - S-i) *| , 

where * = ($"^(2x1), . . . , and is an n x n identity matrix. 

The density in the especial case of a bivariate normal copula has the form 

Cp[Ui,U2) , exp( jexp 



2^ 



0:37-- 2 2(1 

3.5.2 t-copula 

The multivariate Student copula or t-copula with two parameters p and v is also inducted in 
the elliptical family copula, it is defined through the multivariate t-distribution as 

where 
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In the bivariate case, 

r^^-^(-) /^^-^(-) 1 ^ x\ - 2px,x, + xl \ ^ 

The corresponding multivariate t-copula density is given by 

Cp,^K,...,^2) = (^^)n/2r(j^) + — - — J where r = . . . , . 

For ly < 3, the variance does not exist and for u < 5, the fourth moment does not exist. 

As z/ ^ oc, the t-copula Cp^iy{ui, . . . , 1/2) is asymptotically equal to the Gaussian copula Cp{ui, • • 



3.5.3 Advantages of copula in Modelling 

Copulas provide greater flexibility in that they allow us a much wider range of possible dependence 
structures. Imagine we have a set of marginals of a given type (e.g., normal). The classical 
representation only allows us one possible type of dependence structure, a multivariate version of 
the corresponding univariate distribution (e.g., a multivariate normal, if our marginals are normal). 
However, copulas still allow us the same dependence structure if we wish to apply it (i.e., through 
a Gaussian copula), but also allow us a great range of additional dependence structures (e.g., 
through Archimedean copulas). These advantages imply that copulas provide a superior approach 
to the modelling of multivariate statistical problems. 
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Copula And Algebra 



Here we present a relation between copula and matrices operations. For more discussion and 
rigorous proof about this section, we refer to [6J. 

4.1 The *-product Of Copulas 

The notion of * product of copula was introduced by Darsow et al. in the context of Markov 
processes. 

dC dC 

Let C = Ciu, v) be a copula and DiC = ^— , D2C = ^— are two-first order partial derivatives. 

ou ov 

Definition 4.1. Let Ci and C2 be two copulas and define the ^—product of Ci and C2 as: 

^1*^2: [0,1]' [0,1] 

{u,v) ^ (Ci * C2) ^) = / D2Ci{u,t).DiC2{t,v)dt 

Jo 

Theorem 4.2. Given Ci, C2 two copulas, the product Ci * C2 is a copula. 

Proof The two first properties about copulas are obvious since, it suffices to apply the definition 
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of the product of copula. To prove the third property, we compute : 



^(Ci*C2) ([^1,^2] X [vi,V2]) 



{{Ci * C2) {U2, V2) - {Ci * C2) {U2, Vi)} 

{(Ci * C2) K, ^2) - (Ci * C2) K, ^1)} 




/ {D2Ci{m,t).D^C2{t, V2) - D2Ci{m,t).D^C2{t, vi)) dt 




[ [{D,C2{t,V2)-D,C2{t,v,))] [D2Cr{u2,t) - D2Cr{ur,t)]dt 



f [Dr {C2{t,V2) - C2{t,vr))] [D2 - dt, 



since ((^2(^,^2) — ^^2(^,^1)) and {Ci{u2,t) — Ci{ui,t)) are both positive, we deduce that 

^Ci*C2 ([^1,^2] X [^1,^2]) > 0, 



Remark For any copula C, 

• n*(7 = (7*n = n, which means that the copula IT is the null element 

• M*C = C*M = C, meaning that copula M is the identity element 

• the * product is an associative operation, 

• the * product of copula is not commutative (for example W ^ C ^ C ^ W). 

This binary operation could be seen as a continuous analog of matrix multiplication. We have 
also W^W = M. 

4.2 Discrete Copulas and bistochastic matrices 

In this part we give only basic definition, property and proposition of discrete copula. A more 
extensive study can be found in [24j, where Kolesarova et al. have introduced the product of 
discrete copulas and their links to bistochastic matrices. 

Let us denote the grid of the unit square by In := {O, ^, |, . . . , l} . 



meaning that the {Ci *C2)— volume of the rectangle [1^1,1^2] x [^1,^2] is positive. 



□ 
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Definition 4.3. A function Cn^rn '-In^'^m 

it satisfies the following conditions: 

For all i e {0, . . . , n} and j G {0, . . . , m} 

1. C,^^(A,0) =0 = C,,^(0,^), 



[0, 1] is called a discrete copula on 1^ x //" 



while for a// i G {1, . . . , n} anof j G {1, . . . , m} 

5. C^,m ^) + C"^,??! ^) ^ C"^,??! ^) + Cn^m ^) • 

Now focusing on the particular case, where n = m and writing Cn^n Cn- 
Property 4.4. y4 bistochastic matrix is a matrix A {^ij)i<i<n-i<j<n ^^^^ ^^^^ 

n n 

for all k ^ {1^ . . . ^n} ^ aij > and aik = akj = 1. 

k=l k=l 

The relation between a bistochastic matrix and the representation of a discrete copula is stated 
by the proposition. 

Proposition 4.5. For a function Cn : XnX^n — ^ [0, 1] the following statements are equivalent: 

1. Cn is a discrete copula; 

2. there is a bistochastic matrix A {^ij)i<i<n'i<j<n ^^^^ ^^^^ i ^ {0, 1, 2 . . . , n} , 

c;?:=af-,^)=-EE«^"^ (4-1) 

^ / k=\ m=l 

The discrete copula C„ could then be rewritten in the term of n x n matrix, with entries given 
by Pq]): 

Definition 4.6. Let C\ and C\ be discrete copulas both defined on Xn x In and let Ai and A2 
be the nx n bistochastic matrices corresponding respectively to them. Then the discrete copula 
Cn associated to the bistochastic matrix A = A1.A2 is the product ofC^ and C^. This product 
is denoted by 

Cn = Cn^ ^n- 

For some special class of copulas there is a relationship between the ^ product of discrete copulas 
and the *product of copulas in the sense of Darsow. 
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Copula for Modeling and Parameters 

Estimation 



Let us consider, the following statistics problem. Suppose that we have for example the following 
set samples {x\,xl} for i = 1, . . . ,r. And we want to choose the copula which fits them. We 
expose here some of the statistical methods used in that case. First we consider the case where 
we have selected a joint probability density function depending on parameters 9 and we want to 
estimate these parameters. We expose these methods : MLE, IFM and Bayesian approach. Then 
we consider the problem of model selection where we want to select between a few copulas the 
one which fits the best data in the sense to be specified and revise different criteria which have 
been proposed for this. 



5.1 Maximum Likelihood Estimation(MLE) 

The Maximum likelihood estimation (MLE) is a statistical method used for fitting a mathematical 
model to some data. Modeling real world data by estimating maximum likelihood offers a way of 
tuning the free parameters 9 of the model to provide a good fit. Commonly, one assumes that 
the data drawn from a particular distribution are independent, identically distributed (iid) with 
unknown parameters. This considerably simplifies the problem because the likelihood can then 
be written as a product of n univariate probability densities. From (I2.17p . if we suppose that the 
copula, the marginal functions and the densities functions dependent on the parameter 9 ( which 
could be a vector of parameters). 

• The likelihood L{9) is written as : 

T 2 

L{e) = ]lc{FM-,o),Hx2-,oy,o)l[Mxy,0), 

1=1 j 

• and computing the log-likelihood C{9) yields to: 

T T 2 

£ (^) = log c (Fi (4 ; ^) , F2 (4 ; ^) ; ^) + log f,{x)-e). (5.1) 

1=1 1=1 j = l 
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This method estimates 9 by finding the value of 9 that maximises C{9). The value of 9 is then 
performed through the following MLE estimator: 



9 argmax>C(^). 



(5.2) 



5.2 Inference Function for Margin (IFM) 

For purpose of algorithm implementation when the number of parameters is large, the IFM is 
mostly applied. Instead of using MLE to estimate in one step the parameter 9, one can use the 
Inference Function for Margin (IFM) method; for more discussion about this method see [21j . 

Basically it is a two steps method: 

First step: One estimates parameters 9j (often a vector of parameters) for each margin function. 
The likelihood, the log-likelihood and the estimator of 9j are successively 

T 
i=l 

T 



9j = argmax£j(0j) for j = 1, 2. 



(5.3) 



Second step: Estimate the parameter 9c of the joint density function 



T 



c{e„ e,, 02) = log f{xi 4; e,, 4, 0,) 



1=1 



Now the IFM estimator of 9r is 



9c = argmax£(0c5 ^i, ^2)- 



(5.4) 



To look how those methods work; one could assume that the two marginal functions are nor- 
mals with known parameters, and that the bivariate distribution is also Gaussian with another 
known parameter. From those true parameters, one could then generate a Gaussian copula data. 
From both MLE and IFM, compute the estimate parameters and compare to the original true 
parameters. This could be done using a comparison table and also through several scenarios with 
any other bivariate distribution functions. See [37j for a full treatment of these topics, and easy 
implementation for MLE and IFM methods, or to set up a table comparison between MLE and 
IFM for different copulas families. 
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5.3 Bayesian Approach 



As in the MLE approach, if we have a set of samples {x\, x^} for i = 1, . . . , T for which we have 
chosen a parametric copula family c{x \ 9) and a likelihood function 

T 2 

fix I 9) = Y[c{F,{x\ I 9),F,{xl I 9))l[f,{x) \ 9), 

1=1 j=i 

and if we also have some prior knowledge on the unknown parameter 9 in the form of a prior 
probability 7t{9), then the Bayesian approach consists in computing the posterior probability 



f{0 I x) = 



7r(^)/(x I e) n{e)f{x I d) 



n{e)f{x I d) (19 



and then choosing an estimate for from this posterior. The general approach is to choose an 
utility function u{9,0), compute its expected value 



u 



(6) = I u{e,e)f{e\x) 



de 



and choose as a point estimator 



6 = arg min | . 



(5.5) 



Interestingly for different choices of u we find different classical estimators for 9. 
The Expected A Posteriori (EAP) estimation: 

u{9, 9) =\\9-9 IP — y 9eap = J 9f{9 \ x) d9. 

The Maximum A Posteriori(MAP) estimation: 

u{9,9) = S{9 -9) — > 9 MAP = arg max /(^ | x) 

In this second case, we can see easily the link with MLE, because 



9 = arg max /(^ | x) = arg max {£(0) + log7r(0)} . 



(5.6) 



where C{9) is the likelihood given in (15. ip . 



5.4 Choosing The Right Copula 

We have seen that copulas provide greater flexibility in that they allow us to fit any marginals 
we like to different random variables, and these distributions might differ from one variable to 
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another. We might fit a normal distribution to one variable and another distribution to the 
second, and then fit any copula we like across the marginals. But one of the great difficulty is 
to find the right copula since there is a huge number of different copulas proposed. The range 
of parameters for each given copula family is different, then it is difficult to make a comparison 
between them. 

The most used procedure of copula model selection is the one that has larger likelihood. There is 
a comparison between different families of copulas (see [32]), made by computing the Kullback- 
Leibler distance between copulas with densities ci, C2 and the expectation value of ci\ 



-1 rl 

K{ci,C2)=E, 



log 



/ / log [ ci(u^ v) dudv. (5.7) 

Jo Jo \c2{u,v)J 



More precisely for symmetry reasons, the Jeffreys' divergence measure defined by 



^(ci,C2)= / / {ci{u,v) - C2{u,v))\og ( ^^j^^^M dudv, (5.8) 

Jo Jo \C2{U,V)J 



was used, since 



J{ci,C2) = K{ci,C2) + K{c2,Ci). 



Referring to [17] where the Bayesian method to select the most probable copula family among a 
given set is discussed. The prior information of choice is based on Kendall's tau (r) measure of 
association, defined for continuous variables X and Y as 



r(X,y) = 4/ / C{u,v)dC{u,v)-l. (5.9) 
Jo Jo 

As shown by Genest and Mackay [T3], for Archimedean copulas with generator 

r(X,y) = 4 f ^dt + l, (5.10) 

Jo 9^ yt) 

Spearman Rho (p) is another measure which uniquely dependent on the structure of C , but not 
on the behaviour of the margins function: 

p{X,Y) = 12[ [ uvdC{u,v) -3 = 12 [ [ C{u,v) dv du - 3. (5.11) 
Jo Jo Jo Jo 

Kendall's Tau measure and Spearman Rho take advantage of the classical Pearson's Rho (or 
Pearson product-moment correlation) which reflects the degree of linear relationship between two 
random variables X,Y with expected values E{X) = fix,E(Y) = /xy and finite nonzero standard 
deviations (Jx,o-y and defined as : 

and we may also write 

E{XY) - E{X)E{Y) 



r{X,Y) = 



^/E{X^) - E\X) ^E{Y^)- E^Y)' 
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For any increasing functions / and g, we have : 

p{X,Y)=p{f{X),g{Y)), (5.12) 
T{X,Y) = T{f{X),g{Y)) (5.13) 
r{X,Y)^r{f{X),g{Y)) (5.14) 

For example, the bivariate Gaussian copula with correlation 9 has r = |arcsin^. There is 
also another interesting measure in the theory of extreme value copula [21j, helping to compare 
different copula families, so called tail dependence measure as we define below: 

Definition 5.1. Let C{u^ v) l — u — v + C{u^ v) be the joint survival function for two uniform 
(0,1) random variables whose joint distribution function is the copula C. If C is such that 

C(u u) 

lim ^ — Xu exists, then C has an upper tail dependence if Xu ^ (0, 11 and no upper tail 

u^l 1 — U 

Ciu u) 

dependence ifXu 0. Similarly, if lim ^ — = Xl exists, then C has an lower tail dependence 

if Xl G (0, 1] and no lower tail dependence if Xl 0. 



Discussion and technical way to make a good choice of copula could be found also in [7j. Other 
approach of choosing a suitable copula to model dependence based on exponential family is 
proposed recently [23j. 

Definitely « How does one choose a copula ? », this is one question among many other, asked 
by Dr. Mikosch in his paper [28j, when he is raising doubt on the fundamental basis of Copula's 
Theory in comparison with the Theory of Stochastic Processes. 

Prof. Genest answered [llj: 



Model selection is a broad question for which a completely satisfying answer does 
not yet exist, even in the univariate case. The same strategies can be used here as 
in many other modeling exercise, i.e choices can be guided by model properties and 
characterisations diagnostics tools, cross-validation, predictive, accuracy , etc. Given 
that copula modeling is still in a relatively early stage of development, we concur that 
much remains to be done in this regard. For a state-of-art illustration of methodology 
currently available see Genest and Favre [12] . 
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Tomographic Image Reconstruction 



Is it possible to see the interior structure of an object without cutting it open ? The answer is 
yes if we can expose this object to a ray (for example X rays) and to measure its interaction 
with it (for example the X rays radiographies). In this section, we give an introduction to this 
important subject of imaging sciences which helps to solve the previous problem and much more 
arising from different area of science. There is a long list of research domain where tomography 
technique is applied, for example in archaeology, biology, geophysics, oceanography, materials 
science, medical imaging, astrophysics. 

6.1 Tomography technique 

More modern variations of tomography involve gathering projection data from multiple directions 
and feeding the data into a tomographic reconstruction software algorithm processed by a com- 
puter. Different types of signal acquisition can be used in similar calculation algorithms in order 
to create a tomographic image. There are several types of tomography technique associated to 
a specific physical phenomenon. The following table give a list of some methods mostly used: 



Physical phenomenon 


Type of tomography 


X-rays 


X-rays Computed tomography (CT) 


gamma rays 


Single Photon Emission Computed Tomography(SPECT) 


electron-positron annihilation 


Positron Emission Tomography (PET) 


nuclear magnetic resonance 


Magnetic Resonance lmaging( MRI) 


ultrasound 


Medical sonography (ultrasonography) 


electrons 


3D Transmission Electron Microscopy (TEM) 



In our case, we focus on X-rays Computerized Tomography (CT) which used x-ray, through the 
Radon transform. In Medicine for example. Computed tomography (CT) is a diagnostic procedure 
that uses special x-ray equipment to obtain cross-sectional images of the body. The CT computer 
displays these pictures as detailed image of organs, bones, and other tissues. This procedure is 
also called CT scanning, computerized tomography, or computerized axial tomography (CAT). 
One vital application is the diagnostic of cancer; CT is used in several ways: 

• To detect or confirm the presence of a tumor; 

• To provide information about the size and location of the tumor and whether it has spread; 
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• To guide a biopsy (the removal of cells or tissues for examination under a microscope); 

• To help plan radiation therapy or surgery; and 

• To determine whether the cancer is responding to treatment. 

There are also different geometrical view about the object. The representation can be in 2D, 3D 
(static tomography) or AD when adding time parameter, leading to a dynamic tomography {AD 
model representation is used also as a heart model in Radiology). 

6.2 Computerized Tomography and Radon Transform 

The simplest and easiest way to visualise this method is the classical system of parallel projection, 
where the data to be collected as considered to be a series of parallel rays, at position r, across 
a projection at angle 9. This is repeated for various angles. 



Basically the idea of the X-ray CT is to get images of the interior structure of an object by X- 
raying the object from many different directions. X rays go in straight lines inside the body and 
its energy is attenuated more and less depending on the density of the matter in his trajectory. 
Attenuation occurs exponentially in tissue : 




Figure 1 : Parallel beam geometry 




(6.1) 



where / is the attenuation coefficient at position {x,y) along the ray path and 

Lr,6' = {^5 y : r = X cos 9 + ysinO} . 
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So the simplest model relating the log ratio ln(///o) of the observed energy / with emitted energy 
Iq to the spatial spatial distribution / of the body is a line integral: 



Pe{r) = -Ml/Io) 



f{x,y) dl. 



(6.2) 



Therefore generally the total attenuation0 p6i(r) of a ray at position r, on the projection at angle 
9, is given by the line integral (16 .2^ . 



6.2.1 Forward problem 

The forward problem is the one dealing with the mathematical expression of the projections. 



Given find the projections po{r) where ^ G [0, tt 



This is a well-posed problem. According to what we mentioned earlier, in the situation when the 
densities function / is known for all positions (x^y), the projections are given by : 



Radon TZ : Pe{r) = 



/(x, y) 6{r — xcosO — ysinO) dx dy. 



(6.3) 



The relation (16. 3p is known as the Radon Transform (or sinogram) of the 2D object f{x,y), 
consists of X-ray projections along all possible lines in the plane. Each line has a specific direction 
and each direction is uniquely identified by the angle 9. 

5 is the Dirac's delta function defined as : 

5{t) = 0, for t ^ 0, and 
= 1. 

) 

The Dirac's delta function has the fundamental nice property that: 

/(r)5(r-ro)(ir = /(ro). 




Geometric transformation 



If we denote with {x,y) for Cartesian coordinates. And we define (r, 5)51 as the Cartesian coordi- 
nates through rotating {x,y) by an angle Q along the counterclockwise direction (see Figure 1 ). 
The transformation between those two coordinates can be described by the following relations: 



cos Q sin Q 
- sin Q cos Q 



For any function ({x^y), the transformation 

(o{r^ s) = ({r cos 9 — ssinO^rsinO + s cos 9) 
denotes the same function in the coordinates (r, 5)51. 



■"^for continuous value of ^ in < ^ < tt and varying r within < r < 00, we will denote P0{r) by p(r, 0) 
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6.2.2 Reconstruction problem 



Using the result for the relation (16. 3p we give the typical problem of image reconstruction. The 
problem is to estimate a multivariate function f{x^y) from its line integrals po{r), that is called 
an inverse problem and can be formulated in the following ways: 



Given pe{r) for different ^ G [0, tt] find f{x,y). 



It is the same problem, when given the Radon transform (or projections) of a unknown object, 
and trying to find the object. 

The projection-slice theorem tells us that if we had an infinite number of one-dimensional pro- 
jections of an object taken at an infinite number of angles, we could perfectly reconstruct the 
original object. 

One obvious solution is to find the analytic expression of the inverse of the Radon transform. 



1 /*c>o 

Inverse Radon 7?."^ : f{x,y) = — - / / 

2vr Jo Jo 



Mr, 0) 



{r — X cos 9 — y^m 9) 



dr d9. 



(6.4) 



In polar coordinates 



Inverse Radon IZ 



1 



Mr, 9) 



Stt^Vo Jo iCcosie 



-dr d9. 



The book [26j page 108, present the proof given by Radon himself from p3] about the inversion 
of the Radon Transform in and further results. 

However, the inverse of Radon transform proves to be extremely unstable with respect to noisy 
data. In practice, a stabilised and discretized version of the inverse Radon transform is used, 
known as the filtered backprojection algorithm which we described below. 



6.3 Analytical Methods 

6.3.1 Backprojection and Filtered Backprojection 

Let us remind briefly some operators used to express the formula of the backprojection and 
the filtered backprojection. For more detail about this approach, see the book « Principles of 
Computerized Tomographic Imaging » |22] . 

Definition 6.1. The Fourier transform (FT) of a function f{x) G C, x G (— oc,oc) is defined 
as 

roo 
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and the inverse Fourier transform (I FT) is given by 

1 

Conditions for the existence of the Fourier transform are complicated to state in general[2j, but 

/oo 
|/(x)| dx < oo. This re- 
-oo 

quirement can be stated as / G L^(M) meaning that / belongs to the set of all functions having 
a finite L^-Norm. 

/oo 
dx < oo. 
-oo 

More generally, it suffices to show / G for 1 < p < 2. 

Definition 6.2. The Hilbert transform is a linear operator which takes a function, s{t), to 
another function, H{s){t), with the same domain. In signal processing, this operator is used to 
derive the analytic representation of a signal s{t). The exact definition of the Hilbert Transform 
using the Cauchy principal value (denoted here by p. v.) is 

1 s(t) 

H : s{t)=p.v.- / ^^dr. (6.5) 

^ J-oo 

Computationally one can write the Hilbert transform as the convolution: 

H: s{t) = — ^s{t). (6.6) 

Tit 

which by the convolution theorem of Fourier transform^, may be evaluated as the product of the 
transform of s(r) with — isgn(r), where: 

+ 1, ifr>0 
sgn(r) = <( 0, if r = 
-1, ifr<0. 

Therefore to get f{x^y) back, from (16. 3p means finding the inverse Radon transform using the 
backprojection (BP) or the filtered backprojection methods, the above operators are applied. 

The backprojection operator is defined as 

1 r 

Backprojection B: b{x,y) = — / p{xco89 + ysinO.O) d9. 

Instead of dealing directly with the expression (16. 4p . one computes the inverse of the Radon 
transform sequentially using the following steps: 



^The Fourier transform of a convolution is the product of the Fourier transforms 
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Derivative V\ p{r, 6) = 



dp{r, 6) 
dr 



Hilbert Transform H: pir\ 9) = - H ^^1^ dr 

^ Jo {r- rO 

1 

Backprojection B: f{x,y) = — / p{xcos9 + ysinO.O) d9. 

27r Jo 

The model of the backprojection is then defined as: 



f{x,y)^BnVp{r,e). 



(6.7) 



Filtered Backprojection (FBP) 




From the above table where we describe the relationship between V, FT and H, defining: 

/oo 
p(r, 9)6-^'''' dr. 
-OO 

Therefore the properties between V, FT, H yields to 



and 



We have demonstrated the following model of the filtered backprojection (FBP)commonly used 
in X-ray CT, 



FT 




Filter 




IFT 




Backprojection 




> 




— > 






B 



f{x,y) 



Section 6.3. Analytical Methods 



Page 35 



which can be rewritten as: 



f{x,y)=BT-^ \u:\Tp{r,d). 



(6.8) 



Notice that * denotes a 2D convolution, it is also shown that: 

1 



Backprojection B: b{x, y) 



f{x,y). 



In polar coordinates 



Backprojection B: b{^, 0) = — * /(^, 0). 



Bridging the gap between Copula and 

Tomography 



After having reviewed copulas and tomography outlined in the previous chapter, we are now able 
to describe our mathematical approach to the problem of tomographic image reconstruction in 
more detail. We consider only the case of two projections: horizontal 9 = and vertical 6 = 7i/2. 



7.1 Horizontal and vertical projections 




^ = 

Figure 2 : Horizontal and vertical projections 



The particular case where we have only two projections 9 = (horizontal) and 9 = 7i/2 (vertical). 
We substitute in (I6.3p to obtain: 



6 = 0, po{r) = jj f{x,y)5{r - x) dx dy 
^ = |> P7r/2(r) = jJ fix, y)5{r - y) dx dy. 



(7.1) 
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Forward problem: Inverse problem: 

Given f{x,y) find /i(x) and f2{y) Given /i(x) and f2{y) find /(x,?/) 



If now, we denote po = fi and ^7^/2 = /2 the following relation also holds: 



^ = 0, hix)= [ f{x,y)dy 

JR 

^ = ^, f2{y) = / f{x,y)dx. 

^ JR 



(7.2) 



It is important to point out that projections (l7.ip are theoretically the same as marginal distri- 
butions (17 .2^ if we assume the positivity and normalisation. This result in the general case was 
shown by Cramer and Wold in 1936, when they inverted the Radon Transform in the context of 
mathematical statistics(see [4j). 

Now we clearly identity the following situation as an inverse problem typically ill-posed: 



Given fi{x) and f2{y) find f{x,y). 



In fact one of the three conditions for a well-posed problem suggested by the French mathemati- 
cian Jacques Hadamard (existence, uniqueness, stability of the solution) is not satisfied [15]. 

It is clear that given and f2{y), there are infinitely many solutions for f{x,y). 

If we look at the Backprojection(BP) solution and the Filtered Backprojection(FBP) solu- 
tion, for = and 9 = 7r/2, we have 



1 r 

f{x,y) = A cos (9 + 2/ sin 6>, 



e)de 



{fi{x) + f2{y)) . 



This implies that the BP solution, resulting from the trapezoidal rule is: 



f{x,y) = l{Mx) + My)). 



(7.3) 
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And the FBP solution is: 

2 Jo x'-x 2 Jo y'-y 



which can also be implemented in the Fourier domain as 

(7.4) 



fix,y) = \ f I a;i I h{ujr)e^^^^ duj^ + \ f | (^2 | f2Me^'''' du^. 



7.2 Looking at Inverse problem in a different way 

Let us consider the following problem in a different way. Between all the solution which satisfied 
the constraints 

Cl: J f{x,y)dy = Mx) 
C2: j f{x,y)dx^My) 
C3: J J f{x,y)dxdy^l. 
Choose the one which minimise a criterion ^{f) such as: 

^iif) = 11 \fix,y)\' dxdy 

^2{f) = j j -f{x,y)\n{f{x,y)) dx dy. 
Interestingly, if we choose ^i{f) we obtain the Backprojection solution 

f{x,y)^Mx) + h{y) (7.5) 

and if we choose ^2{f) we obtain 

f{x,y) = ^exp(-/i(x)) exp(-/2(y)), (7.6) 

where Z is a constant such that C3 is satisfied. 
If we choose to minimise 

W)^ //(/(...) dxd, 

we obtain 

f{x,y)^h{x)My). (7.7) 
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The above ways to look at the same tomographic image problem, we have discussed is based on 
information entropy. We have used the method of Lagrange multiplier £1 (denoted Cg) to optimise 
the solution under criteria Q^i (denoted the L2 — Norm), Q2 (denoted the Shannon Entropy). 

We have also sat the criterion by modifying the KullBack-Leibler distance in order to obtain 

From those preliminary result, we conclude also that the notion of entropy gives some solution 
we have already found and we will find in the next section about copulaEI. 



7.3 Simulation Results Using Copula 

Those results we propose in this section are a starting point for further research in this area. 

In X-ray CT, if we have a large number of projections uniformly distributed without noise in [0,7r] 
angles, the BP and the FBP are good solutions to the inverse problem. But when we have a few 
number of projections, BP and FBP images felt to be sufficient solutions Figure 17711 

The definition and the notion of copula give us the possibility to propose another way to look at 
X-ray CT method. 

Let first consider the case of two projections. In this case, immediately, we can propose a first use 
which correspond to the copula IT. We call this method Multiplicative Backprojection (MBP). 
This name comes naturally if we compare the two equations (17. 3p and (I7.8p . In practice however, 
we have to normalise each projection in such a way that they can be assimilated to a pdf. 

MBP: 



fix,y) = Mx)My) (7.8) 



Figur d7.3l shows two examples of comparisons between BP and MBP on a few simulated case. 
As we can see with only two projections, there is not any hope to reconstruct a complex shape 
object. We need more projections. 

Now let us have a closer look at expression (I2.17p . in an attempt to find a simple expression of 
f{x,y) and to extend (17. BP . We need one more important property of copula. 



Invariance property 



The dependence captured by a copula is invariant with respect to increasing and continuous 
transformations A of the marginal distributions [34j. This means that the same copula may be 
used for, the joint distribution of (Xi, X2) as (A(Xi), A(X2)) or (InXi, lnX2), and thus whether 
the marginals are expressed in terms of natural units or logarithmic values does not affect the 
copula. 



l£^(/,A,) = Ail7,(/) + A2 j f{x,y) dy) + A3 (/2(2/) - j f{x,y) dy) + A4 (l - ff fi^^v) ^v) 

where i = 1, 2, 3, 4 and j = 1,2,3 then solve ^^^^-^^ = and ^^^^-^^ = to find A^. 

' ' ' ' ' df dXi 

^compare ea. (|7.7p to ea. (|7.8p : and also ea. (|7.5p to ea. (|7.3p 
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Figure 7.1: BP and FBP: using 128 and 2 projections. 



The general case and extension of (17. Sp using the invariance property can be formally written as: 
General MBP: 

" " (7.9) 



f{x,y) = fi{x)f2{y) c{x,y). 



• «Choose» any copula density c{x^y) (one of the list we have given) 

• Normalise each projection in such a way to satisfy pe{T) > and J po{r) dr = 1. 

• For each projection, compute a backprojected image, and just multiply them pointwise, in 
place of adding them up. 



Some few results from the copula-tomography package (see Appendix [Aj) to simulate different 
methods of tomographic images reconstruction (BP, FBP, MBP and much more in future) are 
shown in the next section. 

From equation (l7.9p we have generated then we compute the two marginal functions /i 

and /2. Our reconstruction used a test based on choice of copulas densities c{x^y). 

There is a comparative result from only 5 projections in Figure FTSl between BP and MBP. 
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We clearly observe that, in this particular example of simulation using one Gaussian pdlEI and 
four Gaussians pdf's as original image to be reconstructed; the originate(s) Gaussian(s) are 
reconstructed but having small variance(s). 



^see Figure rrsl a). b) and c) 
^see Figure rOl g). h) and i) 
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Original f{x,y) 



Original fix.y) 



Original f{x,y) 




BP 




BP /(a;,!/) 



MBP f{x,y) 



MBP f{x,y) 



Figure 7.2: BP and MBP on three synthetic examples: using only 2 projections. 
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a)Original f{x,y) 



# 



d)Original f{x,y) 



g)Original f{x,y) 




e)BP f{x,y) 



f)MBP fix^y) 




h)BP f{x,y) 



i)MBP f{x,y) 



Figure 7.3: BP and MBP on three synthetic examples: using 05 projections. 
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Conclusion 



In this report, we have started an other mathematical approach from the theory of copula which 
could be used in tomography. The main contribution of this report is to find a link between the 
notion of copulas in statistics and X-ray CT. For this, first we have presented briefly the bivariate 
copulas and the image reconstruction problem in CT. We have also implemented a Matlab code 
for copula and classical tomography method in order to simulate an image reconstruction and 
then present some preliminary result. 

We could make a link between the two problems of 

i) determining a joint bivariate pdf from its two marginals and 

ii) the CT image reconstruction from only two horizontal and vertical projections, by emphasising 
that in both cases, we have the same inverse problem of determining a bivariate function (an 
image) from the line integrals. 

There are clearly many questions that we have left unanswered, among those questions we draw 
attention to: 

1) the way to develop a strong statistical framework and methods including mixture of copulas 
and taking account of parameters in order to reconstruct with accuracy any complicated shape 
of objects, 

2) how to choose and/or construct a new family of copula having a nice property for image 
reconstruction. 

Further research may be undertaken in order to deepen understanding the relation between copula 
and tomography for applications. Nevertheless, we hope we have managed to give the reader 
some useful insight into, and sparked their interest in this discussion situated in the border of 
statistical mathematics and the fascinating area of imaging sciences. 
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Appendix A 



Tomography-copula package 



Here is a short description of the package we have written using Matlab which allows users to 
simulate a tomographic image reconstruction via a wide range of copula family and mixture of 
Gaussians pdf 's. We hope that this package (still under development) will help to enjoy the joy 
of copula in tomography. 



Copulas Tomography 
Simulate dataft<,y), f_lti) andf_2(y) 
Back- Project ion (BR): fht<Al=f_lt<)+f_2(y) 



Filtered BR (FBP): fh(x,y) = df_l(x)+df_2(y) 



Multiplicative BP (MBR): fht>;,y)=f_l (x).f_2(y) 



Multiplicative FBP(MFBP) 



Mixture of Gaussians fitting 



l^jjj ^ fht<,V)=f_lC^).f_2(y). c(x;y^^ 



IS 

File Edit View insert lools Desktop Window Help 




> 



i 



File Edit View insert lools Desktop Window iHelp 



□ 1^ Q a 




20 40 00 80 100 



i 



File Edit View insert lools Desktop Window IHelp 

D y a , I ^ ^ r? ® I v;- \\u\m\ 

My^.Vl and fhj(i/) -4 




File Edit View insert lools Desktop Window jHelp 

□ y a, ^ ^ r? ®rH5|i_a©i 




File Edit View insert lools Desktop Window IHelp 

DG^yaife|^^n®j45| UM\ 

fhj(K) and fj(x) 




20 40 60 80 100 



File Edit View insert lools Desktop Window iHelp 

DG^ya|fe|^^n® h-g | □ B 

fh,(K) and f,<x) 




Figure A.l: «Menu» showing BP simulation: 4 Gaussians mixtures 
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A.l How it works 
A. 1.1 Menu 

The main menu is titled : Copulas Tomography 

It gives a user interface that allows to make selections and choices from the following preset lists: 

1. Simulate data f{x,y), /i(a:) and /2(y) 

2. Back-Projection (BP): fh{x,y) = h{x) + f2{y) 

3. Filtered BP (FBP): fh{x,y) = dh{x) + df2{y) 

4. Multiplicative BP (MBP): Mx,y) = h{x).h{y) 

5. Multiplicative FBP (MFBP): fh{x,y) = dfi{x).df2{y) 

6. Mixture of Gaussians fitting 
7- fhix,y) ^ fi{x).f2iy).c{x,y) 
8. Quit 



A. 1.2 Menu List 

Let us explain each part of the previous list: 
1. 



Simulate data f{x,y), fi{x) and f2{y) 



In fact, this list is a pop-up menu which contains a group of five choices in its own window. 
Those choices are related to the kind of data to be simulated. Those data represent f{x,y), 
the original image. We offer to the user the choice between one Gaussian or two to nine 
mixture of Gaussians. After the selection of the data, clicking on «quit» will close the 
pop-up menu and allow the user to return to the main menu. 



2. 



Back-Projection (BP): fh{x,y) = h{x) + f2{y) 



This list offers a visualisation of the reconstruction using Back-Projection method. fh{x,y) 
is the image reconstructed. 



3. 



Filtered BP (FBP): U{x,y) = d^jx) + df2{y) 



This list offers a visualisation of the reconstruction using Filtered Back-Projection method. 
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Multiplicative BP (MBP): h{x,y) = h{x).h{y) 



This list offers a visualisation of the tomographic image reconstruction using Multiplicative 
Back-Projection method. 



5. 



Multiplicative FBP (MFBP) Ux,y) = df,{x).df2{y) 



This list offers a visualisation of the reconstruction using Multiplicative Filtered Back- 
Projection method. 



6. 



Mixture of Gaussians fitting 



Source code under development. 



7. 



fh{oc,y) = fi{x).f2{y)Aoc,y) 



Click on this list of the main menu yields to a large number of copula family c{x^y) to 
reconstruct the original object. We offer the reconstruction through Archimedean families 
of copula, discrete copula and much more. One could also add many other families of 
copulas. 



8. Quit : To exit from the main menu. 
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How do copulas look like ? 



Here are some graphical representation of copulas, we have discussed in this report. We remind in 
2D case, that and F2{y) are the marginal cumulative distribution functions (cdf's) related 

to joint cdf F{x,y) and the marginal probability density functions and f2{y) are linked to 

their joint probability density function f{x,y) via the horizontal and vertical line integrals. 
We have also c{x,y) for the copula pdf, to be distinguished from the copula cdf C{x,y) (with 
capital letter C). 






F(x,2/),Fi(x),F2(2/) 



f{x,y) with and f2{y) 



C{x^y) contours plot 




0.2 0.4 0.6 0.8 1 



c(x, y) contours plot 





y) mesh plot 



c(x, y) mesh plot 



Figure B.l: Gumbel copula with a = 1.1. 
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0.2 0.4 0.6 0.8 1 



c(x, y) contours plot 

Figure 




f{x,y) with fi{x) and f2{y) 








C{x^y) contours plot 








c(x, y) mesh plot 



C{x, y) mesh plot 
B.2: Clayton copula with a = —0.1. 




Figure B.3: Ali-Mikhail-Haq copula with a = 



-1. 
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I 




c(x, y) contours plot 




f{x,y) with fi{x) and f2{y) 




C{x, y) mesh plot 




0.2 0.4 0.6 0. 



C{x^y) contours plot 




c(x, y) mesh plot 



Figure B.4: Frank copula with a = 0.1. 




0.2 0.4 0.6 0.8 



F{x,y),Fi{x),F2{y) f{x,y) with and f2{y) C{x,y) contours plot 




c{x^y) contours plot C{x^y) mesh plot 

Figure B.5: Gaussian copula with p = 0.1. 



c(x, y) mesh plot 
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I 



F{x,y),Fi{x),F2{y) 



^1 



0.2 0.4 0.6 0.8 



c(x, y) contours plot 




I 



f{x,y) with fi{x) and f2{y) 




C{x, y) mesh plot 




0.2 0.4 0.6 0.8 1 



C{x^y) contours plot 




c(x, y) mesh plot 



Figure B.6: Farlie-Gumbel-Morgenstern family with a = 0.5 




Figure B.7: Copula with cubic section with a = 1 and /? = 0. 
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